Stochastic method for in-situ damage analysis 
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Based on the physics of stochastic processes we present a new approach for structural health 
monitoring. We show that the new method allows for an in-situ analysis of the elastic features 
of a mechanical structure even for realistic excitations with correlated noise as it appears in real- 
world situations. In particular an experimental set-up of undamaged and damaged beam structures 
was exposed to a noisy excitation under turbulent wind conditions. The method of reconstructing 
stochastic equations from measured data has been extended to realistic noisy excitations like those 
given here. In our analysis the deterministic part is separated from the stochastic dynamics of the 
system and we show that the slope of the deterministic part, which is linked to mechanical features 
of the material, changes sensitively with increasing damage. The results are more significant than 
corresponding changes in eigenfrequencies, as commonly used for structural health monitoring. 
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INTRODUCTION 

It is a crucial task to achieve an early and reliable 
detection of feature-changes of mechanical structures 
caused by damage, fatigue, or other environmental in- 
fluences. Commonly, those detection systems use fast 
Fourier transformation (FFT) to extract system features 
[IH2] and to determine the condition of the system from 
changes in the eigenfrequencies. One drawback of this 
approach is that noisy excitation of the structure broad- 
ens the peaks of the frequency spectrum and thus makes 
it harder to detect changes reliably. In this paper we 
present a method to obtain the dynamical behavior of 
the system and to analyze changes of the system's dy- 
namics due to damages. Our proposed method, which 
is based on the physics of stochastic processes, is suit- 
able for in-situ application, as we show that it is robust 
against changing working conditions. In particular, we 
show how to separate the stochastic response dynamics 
of the system from the deterministic one. Even for differ- 
ent levels of exciting noise or, respectively, of turbulence 
we are able to analyze only the determinism of the system 
dynamics directly linked to the mechanical properties of 
the system. 

The paper is organized as follows. First we present the 
selected experimental system and its numerical model. In 
section three the method and the results are presented. 

EXPERIMENTAL SYSTEM 

As a simple mechanical system we used a one-sided 
fixed beam structure and placed it in the wind tunnel of 
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Figure 1. Schematic of the beam structure. 



the University of Oldenburg (cf . Fig. [T]) . To investigate a 
turbulent flow-structure interaction, a motor-driven gust 
generator was used for the production of the turbulent 
inflow conditions. A sphere was mounted on top of the 
beam to increase the acting drag force. The deflection 
of the structure in the horizontal xy-plane was measured 
with a laser diode mounted in the sphere and aimed onto 
a two-dimensional position sensitive detector (2D-PSD). 
This light pointer principle is known to be highly resolv- 
ing, it is used in atomic force microscopy and has recently 
been used also for new anemometers [4j [5]. Thus, the 
bending of the beam in the range of /im can be resolved 
within /xsec. 

For our measurements two steel beams were prepared. 
While one beam remained undamaged, the dynamical be- 
havior of the other beam was changed in two steps. In the 
first step the beam was treated thermally, it was heated 
at red heat and cooled down fast and in the second step 
the beam was cut at a length of 40 % of its circumference. 
For four different wind speeds (7m/s, lOm/s, 15m/s and 
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Figure 2. Time series of the deflection in y-direction. 



17m/s) time series of ten minutes length were recorded 
with a sampling frequency of 30 kHz. From the measured 
data the deflection of the beam in x- and y-direction was 
calculated. Figure [2] presents a segment of the measured 
time series of the deflection in y-direction. 

In addition to the experiment a numerical model of the 
system has been set up. The beam structure was mod- 
eled by a finite element model in which four elements were 
used, three for the beam and one for the sphere. The me- 
chanical properties were chosen in accordance with the 
parameters of the experimental set-up. To simulate the 
damage the stiffness of the second element was reduced 
in several steps up to 70%. The acting forces were cal- 
culated from a series of Gaussian distributed white noise 
[u x = 7, u y = 13, a x y = 2). The forces were fed into the 
model at the nodes, the distribution was calculated ac- 
cording to the geometry and the drag coefficients. Time 
series of ten minutes length where generated at a sam- 
pling frequency of 30 kHz in order to match the parame- 
ters of the experiments. 



STOCHASTIC ANALYSIS 

A wide range of dynamic systems (in particular if fluc- 
tuating, noisy forces are involved) can be described by 
stochastic differential equations, namely, the Langevin 
equation [5] 



The time derivative of the system variable X(t) (X £ 
]R n ) can be expressed as a sum of a deterministic part 
and the product of a stochastic force T(t) and a 
weight coefficient D^ 2 \ For an ideal process the stochas- 
tic force r(t) is white noise with zero mean, i.e. it is 
(^-correlated and Gaussian distributed. In Eq. (JfJ the 
symbol [-]y refers to the element of the resulting 

matrix. Throughout the paper, we apply Ito's interpre- 
tation of stochastic integrals (cf. [5]). 

To show how our experimental situation can be linked 
to such a stochastic differential equation, we start with 
the idea that the system variable X{t) is the deflection. 



As we are interested in the dynamics we propose that for 
X(t) a general differential equation 



X = f(X(t),u(t)) 



(2) 



holds. Here / denotes an unknown function character- 
izing the dynamics, which depends on the deflection X(t) 
and the wind velocity u(t). The turbulent wind velocity 
acting on the structure can be split up into the sum of 
the mean wind speed u and its fluctuations u'(t) 



u(t) = u + u(t) with (u'(t)) = . 



(3) 



The temporal development of X(t) is obtained by in- 
tegration 



X(t + T)-X(t) 



l+T 



f(X,u(t))dt. (4) 



Using the linearization of /(A, u) with respect to u 



8f(X,u) , df(X,u)dX 



OX du 



f(X,u) = f(X,u) + u'(t) 

OU UA U'U 

(5) 

and that /(A, u), i£ and Jj^ ^ are slowly varying for 
small r, Eq. Q can be written as 



X(t+r)-X(t) = r-f(X,u) + 



df + df_ dx 

du dX du 



l + T 



u'(t) dt. 
(6) 

Analyzing experimental data, the mean value of (X(t+ 
t) — X(t))\ X (t)=x can be estimated by taking all X(t) 
values which are in a close neighborhood to a chosen 
value x. This conditional mean is 



(X(t + T ) - X(t))\ x{t)=x = t ■ f(x, u) , (7) 



X t {t) = D?\x{t)) +J2WD(*)(X(t))] ij T j (t) . (I) using 



(r • f{X,u))\ x(t)=x = t ■ f(x,u) 
and, using {u' {t))\ X (t)= x — according to Eq. d3l, 



df(x,u) df(x,u)dX 



du 



dX du 



t + T 



u'{t)dt 



X(t)=x 



df(x,u) df(x,u) dX 
du dX du 



t + T 



{u'(t))\ x(t)=x d< = 0. 



The connection to the Langevin equation (JIJ can 
be seen by expressing the drift and diffusion terms as 
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Table I. Percental values of the slope of the drift function in in- 
direction conditioned on x — (normalized to the undamaged 
values for each wind speed). 





undamaged 


heated 


cut 


7 m/s 
10 m/s 
15 m/s 
17 m/s 


100.0% ± 1.2% 
100.0% ± 0.4% 
100.0% ± 1.0% 
100.0% ± 1.8% 


94.3% ± 1.4% 
93.7% ± 0.7% 
96.2% ± 1.5% 
97.3% ± 1.9% 


74.1% ± 1.5% 
71.8% ± 1.3% 
71.8% ± 0.8% 
68.6% ± 0.6% 



Kramcrs-Moyal coefficients [5], which are the values of 
the conditional moments for r — > 

DP(x) = Hmi^t + T)-^))^ (8) 

r— >u t 

D§\x) = Bmi((X,(t + r) - x^X^t + r) - Xj ))\ x(t)=> 

Siegert et al. [7] and Friedrich et al. [5] developed 
a method to reconstruct drift Z)W and diffusion 
(Eq. Q) directly from measured data for stationary con- 
tinuous Markov processes (for further details see also [9]). 

Equation should be interpreted in the way that 
for every time t where the system meets an arbitrary 
but fixed point x in phase space, X(t + r) is determined 
by the deterministic function D^{x) and the stochastic 
function ^ D< 2 ) (x)T(t). Both, D^(x) and D i - 2 \x), are 
constant for fixed x. 

To analyze the dynamical behavior of our set-up only 
the deterministic part of the Langevin equation is needed. 
From the derivative of Eq. Q we see that the func- 
tion f(x,u) corresponds to Z^ 1 ' and, more important 
for practical purpose, it is only essential to require that 
the mean of the fluctuations will vanish. This is a much 
weaker requirement than requiring delta correlated and 
Gaussian distributed noise. 

For our measured data we proceed as follows. The x- 
and y-coordinates of the deflection of the beam span the 
phase space of X(t) which was divided into 40 equidis- 
tant bins in each direction (fixing different a;- values). 
The drift function was calculated with Eq. ^ point- 
wise for each bin resulting in two 40 x 40 matrices for 
x- and y-direction. Figure [3^, shows a cut in y-direction 
through the drift function for the main flow direction (y- 
direction) at u = 10 m/s. The slope of the heated beam 
is 6.3% ± 0.7% smaller than of the undamaged one while 
the slope of the cut beam is 28.2%±1.3% smaller. Larger 
errors in the outermost bins are due to the small num- 
ber of events in these bins as large deflections are not as 
frequent as small ones. Table [T] shows the change of the 
slopes of the drift function for all wind speeds. The val- 
ues are presented in percental values normalized to the 
undamaged values for each wind speed. 

From the numerical model simulations a similar be- 
havior was found. Figure [SJd shows that an increasing 



reduction of the stiffness leads to a decrease of the slope 
of the drift. For a reduction of the stiffness by 30% the 
slope decreases by 11.6% ± 0.8%, for a reduction of 50% 
the decrease is 22.9% ± 0.4%. 

The change of the slope can be made plausible when 
one links the drift to mechanical features of the mate- 
rial and interprets it as an indicator how fast the beam 
returns to its position of rest. Here it might be noted 
that the negative slope of corresponds to an at- 

traction to a position of rest defined by the fixed point 
D^(x = 0, y) = 0. A decreasing restoring force then 
results in a decreasing slope of the drift. The cut in the 
beam means a major decrease of its stiffness thus the 
slope of the drift should be significantly smaller. (D^ 
is asymmetric for y < and y > 0, see Fig. [3^,, this is 
likely to be due to the fact that a cut in a beam leads 
to asymmetric weakening of the material, whereas in the 
numerical model the weakening in the volume element 
was symmetric, compare Fig. ^p.) 

As the mean wind speed is not constant in free field 
conditions it is important to know if the slope of the drift 
changes with respect to the wind speed, or more gener- 
ally speaking how depends on u. This is particularly 
important for our case where we only condition on bins in 
x and thus the conditional moments of Eq. ([8| sum over 
different u. The obtained results are shown in Fig. |5J: 
and in Table |T] We conclude that independent of a ques- 
tionable or weak dependence of the slope of on wind 
speed the effect of the damages can be clearly quantified. 

As a next aspect we compare our results with the com- 
mon damage detection method using power density spec- 
tra. Figure [4^, shows the power density spectra of the de- 
flection in y-direction for the different beams. The first 
eigenfrequency of the beam structure is very pronounced 
and the peaks are quite broad due to the noisy excitation. 
The shift in the first eigenfrequency for the heated beam 
is only 0.6 % which is almost undetectable due to the 
broadened peak. The shift in frequency for the cut beam 
can be detected more clearly, it is 13.5%. Compared to 
the relative changes in the slope of of several percent 
for the heated and around 30% for the damaged beam 
the frequency analysis is less sensitive. 

A comparison with power density spectra for different 
wind speeds shows that with increasing wind speed more 
and more eigenfrequencies get excited (cf. Fig. [IJs) . Tak- 
ing also Fig. |4^, into account, one can see that a damage 
excites higher eigenfrequencies in a similar way. (Ad- 
ditionally we realized that higher harmonics are excited 
too, if the measured signal saturates, e.g. by overload- 
ing the sensor or by bounded deflection of the structure.) 
The analysis of higher harmonics of the power spectra to 
quantify damages seems to become quite complicated un- 
der such conditions. From a point of view of system dy- 
namics it is well known that damped relaxation systems 
like a beam structure may perform quite difficult non- 
linear response dynamics under noisy excitations, which 
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Figure 3. (Color online) Drift function in y-direction conditioned on x = (a) for three damages at u — 10 m/s (measurement), 
(b) for three damages (simulation) and (c) of the cut beam for different wind speeds. 
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Figure 4. Power density spectra for the deflection signals in y-direction (a) 
at u = 10 m/s for different damages and (b) of the cut beam for different wind 
speeds. 
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Figure 5. Change of the slope of the 
drift function for increasing damage (nu- 
merical simulation) marked by open cir- 
cles. Corresponding damages of the ex- 
periments marked by open triangles. 



even may become chaotic. Consequently an analysis by a 
power spectrum of such systems is less appropriate than 
getting access to the underlying deterministic part of the 
response dynamics. We want to stress the point that 
our analysis by Kramers-Moyal coefficients is local in the 
stochastic variable x whereas power spectra are global in 
the sense that the full range of the phase space variable 
is processed. 

At last, a first estimation of possible resolution of dam- 
ages by our method is given. From simulations with sys- 
tematic changes of the stiffness the changes of the slope 
of the deterministic part were evaluated as shown in Fig. 
[5| A damage of 5 % is well resolvable. Putting the re- 
sults from our experiments into relation of the numerical 
studies (see triangles in Fig. [5]) the data shows that the 
obtained change in the slope of 6% (cf. Table |l| corre- 
sponds to a (numerical) damage of about 15 %. Thus, 
with our proposed method, it should be possible to de- 
tect even smaller damages of the mechanical structure 



than caused by our heating procedure. 



DISCUSSION 

We showed that analyzing the deflection of a one-sided 
fixed beam structure by means of the reconstruction of 
the stochastic differential equation can be fruitful for 
structural health monitoring. The slope of the drift func- 
tion is a sensitive indicator of the restoring force and thus 
of the mechanical properties. The sensitive detection of 
the drift function enables to show changes in the me- 
chanical material properties and thus enables to detect 
probable damages. Most interestingly the method de- 
pends on noisy excitations caused by the environment. 
Noise helps to enlarge the phase space so that the con- 
ditional moments can be estimated properly. We also 
show evidence that the noise itself will be averaged out 
and will not have an influence on the absolute values of 
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the deterministic drift functions. The ranges over which 
it can be reconstructed will be influenced by the noise 
(cf. Fig. |3jj). In this sense the method is robust against 
changing working conditions. We should point out that 
the deflections of our experimental set-up are in the range 
of /xm, thus a linear response can be assumed. For larger 
deflections further corrections may have to be taken into 
account. Compared to the helpful influence of noise for 
the stochastic analysis, the effect of noise for the common 
analysis by determining the strength of excited eigen- 
modes as peaks in the power spectrum is more compli- 
cated and may even cause less sensitivity. 



A very important application of this method are cases 
where one has no easy access to the considered system. 
We propose that this method can be used for remote 
diagnosis of running, embedded systems. Even chang- 
ing working conditions are likely to be mapped onto the 
noise and thus will be filtered out by analyzing the drift 
function. In principle it should be straightforward to ap- 
ply the method to condition-monitoring of other systems 
and other noisy excitation forces, like for example parts 
of a running machinery. For more complex dynamics like 
chaotic ones a higher dimensional phase space has to be 
used. There are even methods to verify by data analysis 
if such a higher dimensional analysis is required jlO] . 
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